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An efficient numerical method is developed using the matrix product formalism for computing the 
properties at finite energy densities in one-dimensional (ID) many-body localized (MBL) systems. 

Arguing that any efficient (possibly quantum) algorithm can only have a polynomially small energy 
resolution, we propose a (rigorous) polynomial-time (classical) algorithm that outputs a diagonal 
density operator supported on a microcanonical ensemble of an inverse polynomial bandwidth. The 
proof uses no other conditions for MBL but assumes that the effect of any local perturbation (e.g., 
injecting conserved charges) is restricted to a region whose radius grows logarithmically with time. A 
non-optimal version of this algorithm efficiently simulates the quantum phase estimation algorithm 
in ID MBL systems; a heuristic version of the algorithm can be easily coded and used to, e.g., detect 
energy-tuned dynamical quantum phase transitions between MBL phases. We extend the algorithm 
to two and higher spatial dimensions using the projected entangled pair formalism. 

PACS numbers: 71.23.An, 75.10.Pq, 02.70.-c, 03.67.-a 


Classical simulation of quantum many-body systems is 
a major area of research in condensed matter physics [39]. 
In general, quantum systems cannot be efficiently simu¬ 
lated classically for the simple reason that the dimension 
of the Hilbert space grows exponentially with the system 
size. Remarkably, most physical systems are highly non¬ 
generic in that they are local. Indeed, locality is essential 
to the design of efficient classical algorithms. At a basic 
level, any local Hamiltonian with short-range interactions 
satisfies the Lieb-Robinson bound, which states that cor¬ 
relations can only propagate at a constant speed (“linear 
light cone”) [22, 31, 34]. Consequently, the short-time 
dynamics of ID local Hamiltonians allows efficient classi¬ 
cal simulation [38]. The locality of a quantum state can 
be characterized by (i) exponential decay of correlations, 
(ii) area law for entanglement [17], (iii) efficient tensor 
network representation [14, 37, 53], etc. For example, the 
ground states of ID gapped [24] and critical [52] Hamil¬ 
tonians are well approximated by matrix product states 
(MPS) [18, 42] of small bond dimension. As the lead¬ 
ing numerical method for computing the ground states 
of ID local Hamiltonians, the density matrix renormal¬ 
ization group (DMRG) [57, 58] is a variational algorithm 
over MPS. 

MBL is an active area of research studying a class 
of interacting quantum systems that are localized in a 
much stronger sense [1, 35]. In the literature, there are 
several (inequivalent) characterizations of MBL: (I) ab¬ 
sence of transport and vanishing dc conductivity [4, 20], 
(H) logarithmic growth of entanglement with time [3, 26, 
48, 55, 56, 60], (HI) area law for the entanglement of 
and efficient tensor network representation of (almost) 
all eigenstates [6, 12, 19, 47], (IV) (quasi-)local integrals 
of motion [13, 26, 28, 45, 47], etc. One might expect that 
MBL systems allow efficient classical simulation by mak¬ 
ing use of their locality. Indeed, (H) strongly suggests 
that the time-evolving block decimation (TEBD) algo¬ 


rithm [54] can efficiently simulate the long-time dynam¬ 
ics of ID MBL systems, cf. [50] is able to work with up 
to 40 spin-l/2’s. For static properties of MBL systems at 
finite energy densities, exact diagonalization is the only 
prevalent numerical method to date, but it is limited to 
small system sizes [30, 33, 36, 40], and finite-size effects 
are not always negligible, cf. the state-of-the-art exact 
diagonalization is able to work with Hilbert spaces of 
dimension 705432 or < 20 spin-l/2’s [33]. It is highly 
desirable to develop efficient MPS-based algorithms for 
computing the properties of excited states in ID MBL 
systems. 

Ideally, one might wish to compute an exact eigen¬ 
state (up to the truncation of real numbers) rather than 
a superposition or mixture of eigenstates that are close in 
energy. This requires a resolution of the order of the level 
spacing in the energy spectrum, which is exponentially 
small in the system size. We believe that any (possibly 
quantum) algorithm for our purposes satisfies the “time- 
energy uncertainty relation” TAE > I, where T, AE are 
the running time and energy resolution of the algorithm, 
respectively. Thus, a generic (exact) eigenstate cannot 
be efficiently computed even on a quantum computer, 
and cannot be efficiently prepared in experiments if any 
experimental operation allows efficient quantum simula¬ 
tion. In both theory and practice, one could only hope 
to efficiently obtain an approximate eigenstate with an 
inverse polynomial energy resolution. 

Inputing a microcanonical ensemble of inverse polyno¬ 
mial bandwidth, we propose a (rigorous) polynomial-time 
(classical) algorithm that (up to negligible errors) out¬ 
puts a (positive) diagonal density operator (mixed state) 
supported on the ensemble, expressed as a matrix prod¬ 
uct operator (MPO). The algorithm first simulates the 
real-time evolution using the matrix product formalism 
(this step dominates the running time of the algorithm), 
and then suppresses unwanted eigenstates using an en- 
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ergy filtering technique [2, 19, 23]. When acting on an 
eigenstate, the propagator can be viewed as a (clas¬ 
sical) wave whose angular frequency is the eigenenergy. 
In this language, we design an interference pattern that 
selects the eigenstates in an energy (frequency) interval 
of an inverse polynomial bandwidth. 

We present a rigorous and a heuristic version of the al¬ 
gorithm. The former should start with a necessary con¬ 
dition for MBL, i.e., violation of this condition implies 
delocalization. We assume neither (III) nor (IV) because 
these heuristic characterizations of MBL seem to be quite 
strong conditions from a rigorous point of view. Rather, 
we assume that the effect of any local perturbation (e.g., 
injecting conserved charges) is (up to negligible errors) 
restricted to a region whose radius grows (at most) log¬ 
arithmically with time (“logarithmic light cone”). This 
is a minimum condition: Any system with diffusive or 
even ballistic propagation of information can be dubbed 
delocalized. Rewriting this condition as a Lieb-Robinson 
type inequality [19, 21, 29], we simulate the time evolu¬ 
tion efficiently and rigorously using the matrix product 
formalism. In two and higher spatial dimensions, a mod¬ 
ified version of our algorithm outputs a diagonal density 
operator supported on the microcanonical ensemble, ex¬ 
pressed as a projected entangled pair operator (PEPO) 
[51] of quasi-polynomial bond dimension. However, we 
do not know how to compute the physical properties of 
such a PEPO efficiently [46]. 

The heuristic version of our algorithm is slightly dif¬ 
ferent in that it simulates the time evolution by TEBD, 
which is the only possibly non-rigorous step (the energy 
filtering step is always rigorous). Although we do believe 
that TEBD is indeed rigorous (TEBD never gets stuck 
in a local minimum because it is not a local search algo¬ 
rithm), we are not aware of such a proof in the literature. 
This version can be easily coded either from scratch for 
a DMRG engineer or by using well-established numerical 
packages [5, 16]. In practice, it can be used to, e.g., de¬ 
tect energy-tuned dynamical quantum phase transitions 
between MBL phases. Such an unconventional quantum 
critical point was predicted in the random quantum Ising 
chain with weak interactions using the excited-state real- 
space renormalization group (RSRG-X) technique [41]. 
We are trying to observe this quantum critical point nu¬ 
merically (in preparation). 

The quantum phase estimation algorithm is an im¬ 
portant subroutine with various applications in quantum 
computing. It can be used to study MBL [7]. In ID MBL 
systems, it can be efficiently simulated by a non-optimal 
version of our (classical) algorithm. 

Note added .—A recent article [43] appeared while the 
present work was being completed. It proposed a heuris¬ 
tic variational approach to approximate the eigenstates 
of a ID MBL Hamiltonian based on the assumption that 
all eigenstates can be mapped to the computational basis 
states by a single local quantum circuit of small depth. 


An important open problem is to identify the class of ID 
MBL systems in which this assumption is valid. Does it 
break down in systems with localization-protected topo¬ 
logical order [27]? 

A fresh article [59] appeared on 18 Aug 2015 after a 
manuscript of the present work was ready to be but not 
yet submitted to arXiv. It introduced the spectrum bifur¬ 
cation renormalization group (SBRG) as a generalization 
of the real-space renormalization group. Various features 
of ID MBL systems (and the quantum phase transitions 
between them) have been obtained numerically using the 
SBRG technique. We like this technique very much. 

It should be clear that our approach is fundamentally 
different from those of [43, 59]. Eurthermore, the energy 
resolution of our algorithm is 1/ poly n, where n is the 
system size. The energy resolutions of the algorithms in 
[43, 59] appear to be c^/n, where c is a small but non- 
negligible constant. Note that the standard deviation of 
energy for a generic state is &{y/n) due to the quantum 
central limit theorem [11]. 

Algorithm .—The main idea is to annihilate unwanted 
eigenstates by destructive interference. Suppose we are 
working with a chain of n spins, and consider a micro- 
canonical ensemble with an energy interval [E — e, E -\- e\ 
whose bandwidth 2e = 1/polyn is an inverse polynomial 
in the system size. Let P be the projection onto the sub¬ 
space spanned by all states in the ensemble. We choose 
{fj € C, tj G R} for j = 0,1,..., V — 1 such that (i) 
N = polyn; (ii) \tj\ = polyn such that each can 

be computed efficiently; (hi) PQP « Q such that Q is 
almost supported on the microcanonical ensemble; (iv) 

N-l 

Q = E (1) 

i=o 


such that Q is diagonal in the eigenbasis of Et. From 
a technical point of view, we prefer to express Q as an 
integral 


/ + 00 

fit)p(H-E)tdt^ ( 2 ) 
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where / : R —?• C is called the filter function [2, 19, 23]. 
We require (i) lim|(|_>oo f{t) decays rapidly such that the 
integral (2) can be truncated to [t] < T for T = polyn; 
(ii) f*{t) = /(—t) such that Q = is Hermitian. 

Filter functions have been systematically constructed. 
The simplest and most efficient one in the present context 
is a Gaussian filter: 
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Let ||X|| and ||X||i = trVx'^X denote the operator 
norm and the trace norm, respectively. Q is almost sup¬ 
ported on the microcanonical ensemble in the sense that 
\\PQP — Q\\ < Thus, \\PQP — Q\\ is negligible for 

q = polyn. It is easy to see that \\PQP — Q\\i is also 
negligible for q = poly n. The truncation (the second line 
in (3)) error is upper bounded by 

VmJr = (4) 

which is negligible for T = polyn. The discretization 
(the third line in (3)) error is polynomially small for r = 
1 /polyn. 

Heuristically, the propagator can be computed us¬ 
ing the TEBD algorithm [54]. As entanglement grows 
logarithmically with time in ID MBL systems [1, 3, 26, 
48, 55, 56, 60], we expect that is well approximated 
by an MPO of bond dimension = polyt. As¬ 

suming that the running time of TEBD is a polynomial 
in the bond dimension, can be computed in time 
polyt. Overall, our algorithm outputs an approximation 
Q' of Q in time (T/T)polyT = polyn. The correctness 
of our algorithm can be efficiently verified numerically 
by checking whether the energy distribution of Q' is suf¬ 
ficiently sharply peaked. 

If entanglement grows polylogarithmically with time 
(e.g., the critical random quantum Ising chain with weak 
interactions [56]), we expect that is well approxi¬ 
mated by an MPO of bond dimension In such 

systems, our algorithm runs in quasi-polynomial time. 

The operator Q can be viewed in two complementary 
ways. First, Q/ iiQ is a density operator (mixed state) 
almost supported on the microcanonical ensemble. Sec¬ 
ond, Q is an approximate projection onto the ensemble. 
To prepare a superposition (pure state) of eigenstates in 
the ensemble, the Hrst heuristic approach is to maximize 
(f/'IQI^) variationally over all MPS \4’) of small bond di¬ 
mension. The bond dimension should be increased until 
the objective function saturates to 1. The second heuris¬ 
tic approach is to start with a random product state JV'o)) 


which is trivially an MPS. In the jth (j > 1) step, an 
MPS IV'j) is obtained by compressing the MPS 
Repeat these steps until the energy distribution of ]?/;) is 
sufficiently sharply peaked. 

To accelerate our algorithm, we can start with a small q 
such that Q (3) suppresses unwanted eigenstates moder¬ 
ately but not significantly. Then we proceed by repeated 
squaring. In the jth step, we obtain an MPO . We ex¬ 
pect that the bond dimension of grows exponentially 
but not doubly exponentially with j as Q is constructed 
from a MBL Hamiltonian H. 

Proof. —We now simulate the time evolution efficiently 
and rigorously. Let H = Hj, where < 1 acts 

on the spins j and j -f 1. Suppose H is MBL in the sense 
that the effect of any local perturbation is exponentially 
suppressed outside a region whose radius grows (at most) 
logarithmically with time: 

\{'f\A{t)\'iP) - {if\U^ A{t)U\f;)\ < e-f 2 (AAC/)) poiy^, ( 5 ) 

Here, j^) is an arbitrary initial state; A{t) := e“®'^‘Ae®^* 
with jjAjj < 1 is a local observable in the Heisenberg pic¬ 
ture; {7 is a local unitary operator describing the pertur¬ 
bation; d{A, U) is the distance between A and U. As any 
local operator B (with bounded norm) can be expressed 
as a linear combination of local unitary operators, (5) 
simplifies to a Lieb-Robinson type inequality 

II [41(0,5]II poly t. (6) 

Reference [29] proved (without using the matrix product 
formalism) that ( 6 ) implies (i) starting from a product 
state, the entanglement entropy grows at most logarith¬ 
mically with time; (ii) the expectation value of a local 
operator {0{t)) can be computed in time poly(nt/i5) on 
a classical computer, where 5 is the desired precision. 

We now prove a stronger statement: ( 6 ) implies that 
e®'^* is approximated by an MPO of bond dimension 
poly(nt/(j); furthermore, such an MPO can be efficiently 
constructed. Let / be an interval, and define 


Hi = Hi{t) = e-®ffi'=-c<=+i)‘iy^e®'^<'=-'-'=+')‘, 

3&1 


Hl{t) = 


(7) 


( 6 ) implies 

\\Hk{t) - Him = \\UHt)HkU{t) - Hfcll < 


-(m(r)Hfc5(r)) 


dr < / ||[5fc(r),5(o,fc-;] + 5[fc+/_„)]||dr 


( k-l - 1-00 \ ,.t 

j g-f^(li-fel) / poly rdr = poly t, I 7 (t) := 

j=-oo j=k+lj 0 


( 8 ) 


Similarly and hence, 

\\Hiit) - Him < \\Hkit) - Him + \\Hk{t) - Him < polyt -h polyt = polyt. (9) 
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Let T be the time-ordering operator. Following [38], we consider the unitary operators 

Vfc(t) := ^ Vl{t) := Te^-^o 

Vl{t) acts on the spins k — l + \,k — l + 2,...,k + l, and is an approximation to Vk{t). Lemma 1 in [25] and (10) imply 
\\Vk{t) - Vlit)\\ < polyt ^ ( 11 ) 

for I = 0{\og{tn/6)). Iterating this procedure for k = 2lk' with k' = 1,2,..., we obtain 


The rightmost side is a quantum circuit of depth 2, where 
each quantum gate is a unitary operator acting on 0(1) 
contiguous spins. Thus, we have explicitly constructed 
an MPO of bond dimension e*^*^*^ = poly(tn/(5) that ap¬ 
proximates to error 5. 

Classical simulation of quantum phase estimation .— 
Naively, a quantum computer obtains a random eigen¬ 
state by measuring the Hamiltonian. However, this mea¬ 
surement is difficult to perform as it is nonlocal. The 
standard approach is to estimate the phases of [/ := 
where 0 < Hr < 2 tt with t = 1/ polyn (n is the system 
size) prevents the eigenvalues of U from wrapping around 
the unit circle in the complex plane. The propagator U 
can be constructed using the Trotter decomposition [32] 
or more advanced methods [8-10]. 

Let {Ekjl'fk}} be a complete set of eigenvalues and 
eigenstates of H, and suppose we are given a state \ip) = 
(^kl'f’k)- In phase estimation, we add an auxiliary reg¬ 
ister of dimension T/r = poly n (or log 2 (T/t) € Z ancilla 
qubits), and label its computational basis by {\j)}J^Q ^■ 
We (i) prepare the state \ j) normalization factor 
is neglected for simplicity); (ii) apply controlled-H gates 
Tjr times; (iii) apply the inverse quantum Fourier trans¬ 
form |j) —>■ (^’^) measure the an- 

cillas: 

E E ^^1^)1-?') = E E 

i i k j 

r{EkT-2TrJ t/T)\J) 

k J 

^'^Ckr(EkT -2TTjT/T)\'tpk)\J), (13) 

k 

where r(x) := — l)/(e“ — 1) is peaked at a; = 0 

with (i) lima;^o?’(a^) = T/t] (ii) maxc<| 2 ,|<i |r(a;)| ~ 1/c 
for r/T ^ c <C 1; (iii) |r(x)| ~ 1 for generic x. Thus, 
the quantum phase estimation algorithm effectively con¬ 
structs the operator R := r{HT— 2 ttJt/T)t/T for a ran¬ 
dom J = 0, 1,...,T/t — 1. Indeed, it fits the framework 
(1) with a trivial filter function: 

N = T/r, E = 27rJ/r, t, = jr, /, = 1, Vj. (14) 


As such, the quantum phase estimation algorithm can be 
efficiently simulated classically in ID MBL systems. The 
operator R is almost supported on the microcanonical 
ensemble with energy [E — e,E + e] in the sense that 
||Pi?P — i?|| < l/(er), but the trace norm ||Pi?P — P|ji 
cannot be satisfactorily bounded. Thus, a nontrivial filter 
function can improve the performance of the algorithm. 

Applications .—We now demonstrate the applications 
of our algorithm in practical settings. As Q' is an MPO 
of polynomial bond dimension, its physical properties can 
be efficiently computed. 

The eigenstate thermalization hypothesis (ETH) says 
that eigenstates at the same energy density have (approx¬ 
imately) the same local expectation values [15, 44, 49]. It 
is believed that MBL (extended) systems violate (satisfy) 
ETH. We expect that the local reduced density matrix of 
Q/trQ typically deviates from the identity matrix even 
in MBL systems. Although a microcanonical ensemble of 
inverse polynomial bandwidth has an exponential num¬ 
ber of eigenstates, the local reduced density matrices of 
these eigenstates are not completely statistically indepen¬ 
dent (and do not average to the identity matrix) regard¬ 
less of whether ETH holds, because the energy density is 
fixed. 

Spectral functions of local operators, which are in prin¬ 
ciple directly measurable in experiments, can be studied 
with our method. As the energy resolution of our algo¬ 
rithm is an inverse polynomial in the system size, we can 
only obtain coarse-grained spectral functions. This al¬ 
ready goes beyond experiments as the energy resolution 
of any apparatus is hnite. 

Consider the weakly interacting random quantum Ising 
chain 

H = J2 ( 15 ) 

2 G Z 

where Jfs are independent and identically distributed 
(i.i.d.), hiS are i.i.d., and J'fs are i.i.d. random vari¬ 
ables. We assume | J'|’s<C | Ji|, |hi|’s. When Ji,hi, J- < 0, 
RSRG-X predicts a critical energy scale Eg. above (be¬ 
low) which (almost) all eigenstates are in the Hilbert- 
glass (paramagnetic) phase [41]. RSRG-X is (mostly) an 
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analytical approach and involves approximations that are 
not fully controlled. Hence, it is desirable that the predic¬ 
tions of RSRG-X are checked against ab initio numerical 
calculations. 

This energy-tuned dynamical quantum phase transi¬ 
tion is characterized by the spin autocorrelation function. 
Let C{t) := {tp\af {t)a^ {0)\'tp), where If/') is an eigenstate 
of H with energy E. RSRG-X predicts C{\t\ —)■ -boo) > 0 
\i E > Ec and C{\t\ -boo) = Q ii E < E^ [41]. 
Furthermore, dynamical RSRG predicts a power-law de¬ 
cay C{t 1) ~ t- 2 A(£;) jf ^ ^ where A(i?) > 0 
is a monotonically decreasing function of E such that 
A(i? —)• Ec) = 0 [56]. We expect that C{t) and C'{t) := 
ti{af {t)af {0)Q)/ tr Q behave similarly if the bandwidth 
e is not too large. We are trying to compute C'{t) up to 
t = polyn and extract its asymptotic behavior convinc¬ 
ingly (in preparation). 

Higher spatial dimensions. —(3) remains valid in di¬ 
mensions higher than one. It appears that we should use 
the projected entangled pair formalism (a generalization 
of the matrix product formalism to higher dimensions). 
If entanglement grows polylogarithmically with time, a 
heuristic counting argument suggests that and Q are 
well approximated by PEPO of bond dimensions eP°^y 
and respectively. In practice, we do not know 

how to implement TEBD efficiently in higher dimensions. 
In theory, the rigorous version of our algorithm can be 
extended to higher dimensions with a little additional ef¬ 
fort, and it runs in quasi-polynomial time. However, we 
do not know how to compute the physical properties of 
such a PEPO efficiently [46]. 

The author would like to thank Bela Bauer, Joel E. 
Moore, and Remain Vasseur for helpful advices. This 
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stitute for Quantum Information and Matter at the Cal¬ 
ifornia Institute of Technology. 
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